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Abstract. The numerical simulation of wave propagation in semiclassical (high-frequency) 
problems is well known to pose a formidable challenge. In this work, a new phase-space 
approach for the numerical simulation of semiclassical wave propagation, making use of the 
smoothed Wigner Transform (SWT), is proposed. There are numerous works which use the 
Wigner Transform (WT) in the study of a variety of wave propagation problems including 
high-frequency limits for linear, nonlinear and/or random waves. The WT however is well 
known to present significant difficulties in the formulation of numerical schemes. Working 
with concrete examples for the semiclassical linear Schrodinger equation it is seen that the 
SWT approach is indeed significantly faster (in a well-defined sense) to work with than the 
WT and than full numerical solutions of the original equation in the semiclassical regime. 
Comparisons with exact and numerical solutions are used to keep track of numerical errors. 
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1 Introduction 

Phase space methods for wave propagation, such as these presented here, are 
in the junction between time-frequency analysis [Hj [71 [8l [T3l [Til [201 [2T| [22| [27] , 
microfocal analysis [5l [161 [13 [231 [Ml [iQj , PDEs [H [1 [9l [TOl [I2] and physics 
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dHl [IS SOI EZl ES] (at least). Distributions from Cohen's class (i.e. Wigner 
Transforms convolved with various smoothing kernels) have been long well 
known and used in signal processing contexts, as semi-empirical ways to 
'clean up' the time-frequency structure of various chirps, e.g. [13]. Microlo- 
cal measures have evolved in the last twenty years, in close connection with 
PDEs, in fact with a very similar job in mind, e.g. [121 [IZj. The need to ex- 
tend powerful spectral-density-based ideas from stationary to nonstationary 
stochastic problems gave birth to various nonstationary spectral densities 
[31 [TOl [m iini [35j. These three hues of work all converge to the Wigner 
Transform (WT), and explain the recent explosion in the literature for WT 
methods for wave propagation (see e.g. [SlEHSElEliiniiniinillEISSlISaESl 
ESI EHl EH ESI ESI SO] ) • One soon discovers however that, despite the analytic 
and asymptotic power of these methods, there has been much less work in 
computing with WTs (for numerical analysis with WTs see [H El II])- This 
is not an accident, and many people complain that it is not at all practical 
to compute with WTs, see e.g. OUT]. Computations with Wigner measures 
are feasible, but this is a different thing, and we discuss it below. In this 
work we propose a way to go around the WTs 'bad' features, and exploit its 
'good' features in the numerical solution of semiclassical wave propagation 
problems. To do that, we work with a particular distribution from Cohen's 
class, which we will call the Smoothed Wigner Transform (SWT). 

This paper is a continuation of [2j. The two papers are complementary 
in the following way: in [2] the results corresponding to individual numeri- 
cal experiments were presented in some detail. The objective there was to 
present direct evidense that 'the SWT approach works'. However, as is dis- 
cussed in section 14. 2[ the efficiency of the method can only be compared to 
the more conventional approaches when working with a population of nu- 
merical experiments for different values of the semiclassical parameter e. In 
this work we focus on this synoptic question, and present clear evidence that 
the SWT approach is indeed more efficient - in a well defined sense - than 
more conventional approaches. 

Historically, the WT 

W^[f,g]{x,k)= I e-'-^'^y f [x + '-^) -g {x - '-^) dy 

was introduced for the study of the semiclassical limit of quantum mechanics 



2 



i.e. the high-frequency Schrodinger equation. Consider the system of 
hnear equations 



£§^u^ (x, t) + U {x, edx) (x, t) = 0, 

for the complex-valued (i-component wavefield {x,t) = {uf {x,t))-^^ ^, 
where x G R" and the d x d matrix-valued Weyl symbol L*^ (x, k) of the 
spatial pseudodifferential operator L"^ {x,edx) is ant i- Hermit ian, L^{x,k) = 

—L^ (x, k) . This is a natural way to quantify the requirement that equa- 
tion ([1]) is a wave equation; indeed many interesting linear wave equations 
have anti-Hermitian Weyl symbols. The use of the Weyl pseudodifferential 
calculus with WTs is natural, see e.g. [HI [33]. The small parameter e scales 
the typical wavelength of the wavefunctions; the length scale at which the 
coefficients vary (or equivalently one over the slopes of the coefficients), as 
well as the typical propagation distance we're interested in are significantly 
larger than e. Many semiclassical problems in quantum mechanics and long 
distance propagation problems in (classical) continuum mechanics can be 
cast in the form ([T]) [T71 [35l |3H] . For the remainder of this paper when we 
refer to a high-frequency wave propagation problem we will mean a problem 
that can be cast in the form ([1]) with anti-Hermitian Weyl symbol. 

There can be (at least) two problems associated with the IVP ([1]): The 
high-frequency-limit problem: recover the limits of certain observables of the 
wavefield, such as the L^-norm density. 



(x,t) := lim 



^< (x,t) 



and the energy density 



(x, t) := lira (L^ (x, edx) (x, t) , tt^ (x, t))^^ . 

On the other hand, the high-frequency problem consists in solving for the 
observables at a given nonzero e <^ 1; naturally, simplifications making use 
of the fact that £ is small can be made, but the idea is not to lose all the 
^-dependent information. The difference is nontrivial; for example it is often 
the case that limit-observables have singularities (become infinite) while for 
any e > the same observables stay finite for all spacetime points. 
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WT methods for high-frequency wave propagation problems consist in 
studying the Wigner matrix {x, k,t) = [W^ [^i^'^'j] (^) ^)] j ^-i ^ i^" 
stead of the wavefields {uf {x,t))^^^ ^ themselves. It can be said that this 
is a good idea for two reasons: first of all, due to 'time-frequency density' 
character of the WT wave-like functions are in principle 'nicely' represented 
in phase space. Moreover, the operator for the evolution in phase space is 
(asymptotically) simpler, and becomes a Hamiltonian flow in many cases. 
In fact this is a general theme in WT-based methods; there are 'two trans- 
forms': one of the 'data' {x,t) {x,k,t), and one of the 'dynamics' 
under which the data evolve. 

A formidable tool for high-frequency-limit problems consists in Wigner 
measures (WM), 

[T6l [26] . Pretty much all bilinear observables of the wavefields can be ex- 
pressed (in the limit) in terms of the WM. Recently computational methods 
for the high-frequency-limit wave propagation based on the WM have been 
developed [231 [Ml [25]. 

Although often discussed in the context of the WT, the question of 
smoothing phase-space dynamic models and its implications has only rarely 
been actually studied in a deterministic context - see however P, [3l] . 

The remainder of this paper is structured as follows: in section 2 we 
review the WM approach for the limit problem; SWTs and the phase-space 
reformulations of problem ([1]) in terms of them are reviewed in section 3. 
Numerical results on the SWT methods for the high-frequency problem are 
presented in section 4; conclusions and further work are discussed in section 
5. 

2 The WM approach for the high-frequency- 
limit problem 

For an excellent (defining) overview of the WM approach one should refer 
to [T7] and for a comparison with WKB-based asymptotic techniques to 
|36j . WMs are a refinement of microlocal concepts related to microlocal 
defect measures, H-measures etc [HI [T71 [26]. They provide a description 
for 'asymptotic signals'; properly speaking, an asymptotic signal is a family 
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of functions {W^(x)}^^q that in some sense share the same time- frequency 
structure. A prototype asymptotic signal is given by the WKB ansatz 

J e>0 

and the corresponding WM is 

(x, k) = \A {x)f 6{k- V^S (x)) , 

unfolding amplitude and 'local wavenumber' (often referred to as 'instan- 
taneous frequency' in time-frequency analysis) information in phase space. 
Naturally, a unique 'local wavenumber' doesn't always exist, as e.g. for 

multi-WKB signals -u^ (x) = ^ v4„ (x) e^'^"*^^\ In fact an aspect of the 

n=l 

WKB ansatz is that the number of 'branches' at each point (x, t) is a 'hid- 
den' parameter. This becomes very important when seeking a solution of ([I]) 

in the form u'^ (x,t) = ^ An {x,t) e^^"^^'^\ because at different points of 

n=l 

spacetime (x, t) there are different unknown numbers of (nonzero) branches 
= N{x,t). The fact that WMs treat this (i.e. the 'number of different 

local wavenumbers') in a more subtle way is an important advantage in the 

study of high-frequency limits [36j. 

For multicomponent wavefields the WM is a positive-definite-Hermitian- 

matrix- valued measure on phase space [17] . To link matrix- valued WMs with 

physical quantities in general one must work in the context of a wave problem 

like ([T]). In that case 

A° (x) = j ^)) dk 

is the total (i.e. taking into account all components of the wavefield) limit 
norm density and 

(x) = j (L° (x, k) (x, k)) dk 



the total limit energy density [17,35]. Observe that the WT allows the con- 
struction of phase space densities, 'elaborating' on the more traditional phys- 
ical space densities, i.e. 

(x, A;) = tr [W^ (x, k)) , (x, k) = tr (L° (x, A;) W'^ (x, k)) . 



The derivation of exact equations for the evolution in time of the WT 
(x, k, t) corresponding to a problem of the form ([T]) is a well studied topic 
[m l26l I2HI [291 EHl EH [39l HO]. From them, asymptotic equations governing 
the evolution of WMs (x, k, t) can be obtained. In simple cases they con- 
sist in decoupled Liouville equations in phase space for scalar phase-space 
measures [T71 Therefore, WT-based numerical methods for the high- 
frequency-limit problem deal with the numerical solution of the Liouville 
equation in phase space when the initial condition is an appropriate mea- 
sure. Moment and level set methods have been formulated and used in this 
connection P [231 [2ll [25] . 

A natural question arises here: can't these methods be used for the small 
5 > problem? So far the answer seems to be negative. The e — > weak 
limit 'kills' a lot of fine scale features. If they are not killed, but included 
in the numerical problem, we get an overly heavy problem, that typically is 
slower to solve than a full numerical solution of ([T|), thus beating its purpose. 
The need to put the extra data under control clearly emerges; this is possible 
in terms of the SWT. 

3 The SWT approach for the high-frequency 
problem 

3.1 The WT vs. the SWT 

The Wigner Transform is defined as a sesquilinear mapping, 

W : f{x),g{x)^W' [f,g] {x,k) . 

When / = (yf it is called the Wigner distribution of /, and denoted as 
[f] {x, k). The Wigner distribution of a wavefunction / (x) is a good 
way to realize a joint breakdown of the wavef unction's L^-norm over space 
X and wavenumber k - with the caution that it takes on negative values as 
well. The books [121 US] are two very important sources on the WT and its 
properties, [13] emphasizing the signal processing point of view, and [IS] the 
relation between the WT and pseudodifferential operators. 

Cohen's class of distributions [3 [H [131 [211 [22] is defined as the class of 
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all sesquilinear transforms of the form 

y,u,z 

for any distribution (pc {z, y) (restricting /, g to test functions). An equivalent 
definition is the class of transforms that results from convolving the WT with 
a distributional kernel 

Kc (x, k) = i^c (2, y)] 

C [/, g] {x, k)= j Kc{x- x', k-k') W [/, g] (x', k') dx'dk'. 

Among them there are many attractive alternatives to the WT. When the 
kernel is chosen to be the Wigner distribution of some function h, Kq {x, k) = 
W [h] {x,k), and if f = g, we get the spectrogram of / with window h. 
Smoothing WTs with an appropriate kernel tames the interference terms, but 
doesn't necessarily kill them completely; there is a balance between smooth- 
ing enough and not smoothing too much. In this connection we define the 
scaled SWT as the sesquilinear transform 

= ^ / e [/, g] {x', k') dx'dk'. ^ ' 



This is a (scaled) WT convolved with a tensor-product Gaussian with 
space-domain variance proportional to eal and wavenumber-domain vari- 
ance proportional to eaf,. The scaling is selected with the problem ([1]) 
and the semiclassical regime in mind; see Figures 1,2 for intuition on the 
scale of smoothing. The SWT provides coarse-sca 1^1 description of the time- 
frequency structure of the underlying signal / wavefield; its resolution, i.e. 
how much information is effectively discarded, can be tuned through the 
smoothing parameters. 

^Different authors favor different terms, e.g. slow-scale, coarse-scale, macroscopic, ho- 
mogenized etc. Each of these terms stands for a description which contains information for 
the underlying object, understood to be sufficient in terms of certain qualitative criteria, 
but which is not complete - usually it is essentially less than complete, thus leading to 
compression as well. 



7 




0.2 0.4 0.6 0.8 



Figure 1: Wigner Transform of the 'tanh chirp', f^{x) = A(x)e^'^'^ A{x) = 
g-25{x-0.5)2^ _ ^l^g (gl0(a-0.5) ^ g-l0(x-0.5)j (Remark: Image quality has to 

do with - the very understanable - arXiv size Hmits, and does not represent the quahty 
of the output; see e.g. 0). 
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Figure 2: Smoothed Wigner Transform of the 'tanh chirp', correspondig 
to the scahng ax^CTk = 0(1) (Remark: Image quality has to do with - the very 
understanable - arXiv size limits, and does not represent the quality of the output; see 
e.g. ID). 
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3.2 Equations for SWTs 

In [2] (an early version of which can be found in the arXiv) we presented the 
derivation of exact equations for the SWT of a wavefunction. That result 
essentially stated that, given the multicomponent problem ([1]), its smoothed 
Wigner matrix satisfies 



e^W + A + B^ = 0, 
at 



where 



(3) 



A = L' 



X - —dk - ^d,, k + -^d, - ] w, 

Am An Am Att 



\ Am An Am An 

with initial data 



i,j=l,...,d 



We have used the abbreviation W 



W-'""'^ [ul,u^] {x,k,t) 



i,j=l,...,d 



and, as we discussed earher, the Weyl PDO calculus. For concreteness we 
will now focus on the (scalar) linear Schrodinger equation on M, i.e. 



{x, t) - i^Au'' {x, t) + iV (x) (x, t) = 0, 



■u^ {x, 0) = Uq (x) 



(4) 



In that case (131) becomes 



f^W{x,k,t)+(2nk^^ 



+ 



2 dxdk 



W{x,k,t) + 



+ lReUv(x- ^d, -'-tdJW (x, k,t))= 0, 



(5) 



W (x, k, t) = W"-^"^'^ [ul] (x, k, t) . 

This is the exact equation, not an asymptotic approximation. Observe 
that the smoothing parameters ax, (Tk are now part of the data of the phase- 
space equation. 'Little' smoothing, i.e. a regularized WT, corresponds to 
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CxiO'k ^ 1 in the scaling ([2]). In that case the asymptotic phase space equation 
is 



In particular, the leading order part is a Liouville equation, and is the same 
as the corresponding leading-order equation for the WT. 

The numerical solution of these equations for the SWT are expected to 
provide a way for a very efficient coarse-scale solver for wave propagation 
problems. Indeed, in the more standard approach of solving the original 
wave equation - in this case equation (jl]) - we need to keep track of the 
wavefunction 'exactly'. That is, even if a coarse-scale amplitude is all we are 
looking for, in order to be able to predict how the coarse-scale amplitude 
will evolve, we need as well to keep track of the fine-scale features of the 
wavefunction - up to numerical errors of course. The SWT approach on the 
other hand, enables us to formulate (exactly or asymptoticly) a numerical 
problem for a coarse-scale description of the wavefield. Thus, we have to keep 
track of significantly less information, and ultimately we are able to carry 
out the computations significantly faster. 

4 Numerical results 

4.1 The numerical scheme 

The subject of this paper is the formulation of a numerical method suitable 
for the simulation of semiclassical problems in terms of the SWT, and the 
investigation of its efficiency as compared to more conventional methods 
for the same problem. The test cases we work with is the I VP (jlj), up 
to time tjnax, independent of e, with initial data scaled with e, for different 
initial conditions and potentials. The model we discretize is the Liouville 
equation for the SWT with Ua;, cr^ ^ 1, i.e. the leading order part of The 
conclusions of this investigation, based on the numerical results of sections 

4.2 and 4.3 are discussed in section 5. An outline of the algorithm follows. 
The WT is computed on a Cartesian grid in phase space with the FFT, 

with complexity O(A^^logA^) |27j. The complexity for the computation of 
the SWT is O [LP'N'^ log A^) where L is the number of sampling points needed 




W{x,k,t) = 0{e^). 



(6) 
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for the smoothing kernel. When a^^a^ = 0(1), L is of 0(1). Adaptive 
computation of the SWT that doesn't spend much time on regions of phase 
space with no energy is also possible |32], and might be essential for two- and 
higher dimensional problems. 

The Liouville equation is solved numerically with the use of particles, i.e. 
the numerical implementation of the method of characteristics (related, but 
not identical, to [H [HI El] ) • An initial population of particles is created, so 
that the SWT can be interpolated up to an error tolerance from its values 
on them. The trajectory of each particle is computed according to Hamil- 
ton's ODEs with a Runge-Kutta solver; the value of the density on each 
particle remains unchanged in time. The solution at each moment in time is 
constructed by interpolating the density from its values on the particles. 

4.2 SWT computations vs. solutions of the original 
problem ([4]) 

The results of the SWT method are compared with exact as well as numer- 
ical solutions of the corresponding Schrodinger equation Being able to 
compute faster with the SWT than by numerical methods for the original 
problem (jl]) (or ([T]) in general), is essential to our point of view. So far, the 
method we present here appears to be significantly faster than WT methods 
and numerical solutions of (jlj); however first we have to make precise what 
exactly we mean by that. 

It is well known that as the semiclassical parameter e becomes smaller, 
the numerical solution of (jl]) eventually becomes intractable. Saying that a 
numerical method is slower or faster than another in the e <C 1 regime we 
mean the computation time for the 'same' problem grows slower or faster 
in e for the one method than for the other. Hence we're looking for the 
semiclassical complexity T = T (e), that is the asymptotic dependence of the 
total computation time on the semiclassical parameter. In the same way, 
we're intersted in the rate of growth of the degrees of freedom needed for 
each numerical method D = D{e). That way we can compare numerical 
methods that are implemented in very different ways, and still get a well 
defined, meaningful result. 

As was shown in ||4j, meaningful numerical solutions by means of finite 
differences methods for can be obtained with two different meshing strate- 
gies: a finer meshing strategy is needed for the correct approximation of 
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the wavefunction, while a coarser meshing strategy can lead to correct slow- 
scale observables (but wrong fast-scale information due to numerical disper- 
sion). As we discussed earlier, the SWT method can only recover slow-scale 
information. Independent numerical experiments, carried out by Kostas Poli- 
tis (NTUA), lead to findings consistent with those of [Ij: the semiclassical 
complexity of standard finite differences methods appears to be O {s~^) or 
more, depending on the meshing strategy. 

A wavelet method developed by K. Politis shows a rate of growth of 
computation time around 0(iV^'^), where N is the number of points (i.e. 
translated versions of the finest-scale wavelets) used; is proportional to £~^. 
Thus wavelet methods seem very promising for semiclassical computations 
as well. Wavelet based solutions also appear to be better behaved than finite 
differences in terms of numerical dispersion, and can be accelerated with 
thresholding. 

In the numerical experiments carried out so far, the SWT approach suc- 
cessfully captures the slow-scale observables of the wavefunction with semi- 
classical complexity no more than O {—e~'^log {e))\ see Figures 3, 4. It is in 
that sense that our preliminary investigation shows that the SWT method is 
/aster than full solutions of the problem (jlj). (In all examples considered the 
SWT was significantly faster than the wavelet method, although the rates 
of growth of computation time in are relatively close for the two meth- 
ods). A more thorough study of this question of semiclassical complexity, 
including more theoretical aspects as well as different numerical methods 
(e.g. time-splitting spectral methods, which seem to be particularly well 
adapted to semiclassical computations) is in order, and in fact in progress, 
in collaboration with other groups as well. 

4.3 SWT vs. WT computations 

WT methods are slower than either full solutions of (jlj) or SWT solutions. 
As compared with the SWT method we propose, the only difference would be 
another, more complicated initial condition. Apart from being represented 
with more points, the initial WT is qualitatively different than the SWT. It 
features many short waves with important amplitude, being a high-frequency 
2-d wavefunction itself. 

Our numerical experiments indicate that WT methods are in any case one 
order of magnitude in slower, and sometimes even more. One particular 
phenomenon we observed in those numerical experiments is that, especially 
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0( epsilon ^ ) 
0( epsilon ^ ) 
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Figure 3: Total computation time against e ^. All problems are spe- 
cial cases of (4); Problem 1 corresponds to u\}^x) 

i {tank {Q.87{x + 2.42)) + 1) {tank (6.87(2.42 - x)) + 1), S{x) = - x^ + 
2x; V{x) = 0; Problem 2 corresponds to Uq{x) = A{x)e^^^^\ A{x) as 
above, S{x) = — ^ + 2x; V{x) = x; Problem 3 corresponds to Uq{x) = 
e~T^^'^ + e~°'ii^' + e~°io/'^^; y(z) = Problem 4 corresponds to 
„^(x) = e-(i+f)^'-2^-4 + e-(i+f)^'-^-i + e-(i+i)^'-i^-t; ^(x) = 0. The 
SWT solutions are found to be valid slow-scale representations of full nu- 
merical solutions of Problems 1,2,3 (courtesy of K. Politis), and of the exact 
solutions of Problem 4. Quantifying in a satisfactory way how good a 'slow- 
scale representation' is has proved somewhat subtle, and is in progress; see 
[2] for qualitative comparisons. The time for the computation of the SWT 
of the initial condition, although asymptotically important, is very small in 
this examples. Certain normalizations were necessary to fit all the data in 
the same axes; in any case the slopes are not affected. 
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Figure 4: Number of particles used against e 



in problems with strong interference terms, in order to recover a meaningful 
approximation of the evolution of the WT many more points are required 
than in order to represent it initially. That is, if we use just enough particles 
to approximate the initial WT well, then often our numerical solution evolves 
to something very noisy; we need to use many more particles than that. This 
point is in contrast with the respective behavior for the SWT, and in fact we 
have no precise understanding at this point of how the number of particles for 
a given WT should be selected. (When this behavior kicks in, the complexity 
of the WT solution seems to be more than O {e'~^)). 



5 Conclusions — Further work 

Semiclassical computations are well known to be a very demanding numerical 
problem. However, it seems that the more precise study of how different 
methods, in general solving for different things, behave as the semiclassical 
parameter goes to zero has started only recently. In this work we carry out 
a preliminary numerical investigation of a recently formulated SWT-based 
method for semiclassical simulations [2], and find that it is in general faster 
than full solutions of the original (physical space) PDE, based on numerical 
experiments, consistently with analysis of the more traditional methods [1] . 
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Also, that - as expected - SWT based computations are significantly faster 
and better behaved than WT computations. 

As we mentioned earlier, equation ([5]) is in fact a whole family of phase- 
space models, which can then be asymptotically treated to create a hierarchy 
of approximations for each phase-space model. In this paper we have only 
worked in a specific asymptotic regime, the semiclassical regime. One reason 
for this is that it allows direct comparisons with existing WT-based stud- 
ies; also the history and evolution of phase-space methods draws heavily on 
semiclassical problems. However it should be clear that the principal rea- 
son for the SWT being 'faster' is the fact that it is, and remains as time 
passes, a much simpler function than the WT. In other words homogeniza- 
tion here comes from the smoothing, and it is an essentially different thing 
than asymptotics. 

These insights lead to the following conclusion: numerical methods cor- 
responding to the exact equations should be formulated and compared to 
full solutions for different problems (e.g. 'any' linear problem) and various 
smoothing choices; we can then look for practical, possibly problem specific, 
'slow-scale solvers'. Thus, further work at hand includes the treatment of dif- 
ferent models included in ([5]). Numerical methods for such models are needed; 
parallel implementation might be crucial for higher- dimensional problems. 
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